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1. Introduction 

When computing the electromagnetic interaction with large complex geometries there is 
generally a trade off between the solution accuracy and the degrees-of freedom representing the 
discrete surface current densities. Most popular electromagnetic scattering codes based on 
method of moment (MoM) solutions employ low-order basis functions and testing functions. 
That is functions that lead to fixed error convergence rates that are typically 0(h 2 ) or less. 

Unfortunately, such techniques are deficient in that: 

• Typical simulations employing 10 to 20 unknowns per linear wavelength provide 3 - 5 % 
accuracy. 

• If increased accuracy is required from the simulation, a significant increase in the 
computational resources is necessary as characterized by low-order convergence. In fact, 
codes employing the popular Rao-Wilton-Gliss (RWG) basis [1] realize only linear 
convergence. That is, a factor of two increase in the number of unknowns only reduces 
the error by a factor of two. 

• Efficient predictions of solution error is not possible. 

In contrast to a low-order scheme, a high-order method is defined as being not just a scheme that 
can provide high orders of convergence, but also a scheme that allows one to control the rate of 
convergence. That is, for a fixed order p, one can decrease the mesh size h realizing a geometric 
rate of convergence 0{hP). Or, one can increase the order p to realize an exponential rate of 
convergence. Furthermore, a high-order technique should be capable of substantially reducing 
the solution error at the expense of only a modest increase in computational resources and 
complexity. 

The locally corrected Nystrom (LCN) method has demonstrated high-order error convergence 
for RCS computations [2, 3] and has been shown to be computationally efficient. The LCN 
method has the distinct advantages that: 

• The error convergence can be increased exponentially by simply increasing the 
quadrature order of the underlying Nystrom scheme, avoiding remeshing for higher 
accuracies. 

• By controlling error convergence, mathematical error estimates can be obtained 
efficiently. 

• The computation of matrix elements is computationally efficient. Typically, this is 
dominated by constructing a sparse local correction matrix, which is 0(N)). 

• Improving solution error typically requires one to simply increase the quadrature order on 
the patches - thus avoiding remeshing. 

• Curvilinear meshing allows for cells that can have linear dimensions that are in excess of 
a wavelength on a side. This is amenable to wide band simulations that allow for a 
coarse meshing that can be applicable to a very broad spectrum of frequencies. Accuracy 
is then controlled by simply manipulating the number of quadrature points. 

• The point based discretization is amenable to fast methods 
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High-order MoM solutions are also realizable by employing higher-order basis and testing 
functions [4, 5]. However, experience has shown that the LCN technique is more efficient and 
of lower complexity than the MoM. The fill time for the MoM solution is 0(N 2 ) and is burdened 
by the double numerical integrations required to compute the entries of the impedance matrix to 
sufficient accuracy. This is in contrast to the single integrations required for the local corrections 
of the LCN method. Furthermore, when exploiting large patches (possibly one or two 
wavelengths along linear dimension), the point-based discretization of the LCN simplifies the 
grouping required by fast methods. 

The focus of this paper is to provide an in depth development of the locally-corrected Nystrom 
scheme as applied to the electric field and magnetic field integral equations. Section 2 outlines 
the general locally corrected Nystrom scheme. Section 3 provides an in depth treatment of 
performing the exact integrations of the kernel convolved with the basis functions for the electric 
field integral equation. Section 4 then treats the magnetic field integral equation. Section 5 then 
details the computational treatment of the local corrections for general basis. Finally, Section 6 
provides some numerical validation. 
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2. The Locally Corrected Nystrom Method 


The Nystrom method [2, 6] is a simple technique that provides an efficient high-order solution 

of integral equations. Consider the integral equation: 

4 >(r) = \G(r-r')J(r')ds' (2- 1 ) 

n 

where J(f') is an unknown current density, 4(r) is the known forcing function evaluated at 
position r , and G is a non-singular kernel. It is assumed that a numerical quadrature rule can be 
defined over the domain Q that can compute the convolutional integral to a desired accuracy. 
Defining the abscissas and weights of a N- th order quadrature rule to be C 2 - 1 ) is 

approximated as: 

< 2 - 2 > 

11=1 

Then, sampling (2.2) at the N abscissa points leads to a square linear system of equations with 
the m-th equation expressed as: 

*(?.) = < 2 - 3 > 

»= 1 

This provides for a solution for J(f) at the discrete points r n '. Assuming that the geometry is 
smooth and G is non-singular, the order of the accuracy would be commensurate with the 
underlying quadrature rule. Thus, if the quadrature rule is exponentially convergent, the solution 
for the unknown current will be as well. 

Unfortunately, the Nystrom method cannot be applied directly to singular kernels. This is true 
for the obvious reasons: I) numerical quadrature integrates singular functions to low order, and 
2) if G(r m - r 7 [) is infinite when r m = r' n , then (2.3) is ill posed. To circumvent this, a modified 

technique known as the locally corrected Nystrom (LCN) method [2, 6] has been introduced. To 
this end, let 


Gw, = 


G(r m ~Z), 

0 , 


m n 


m = n 


be the discrete kernel represented in (2.3). Next, define an exact kernel to be G m „. Then, a local 
correction matrix L m n is introduced such that 

G =G +L . (2.5) 

^rn,n rn,n v 7 

That is, the superposition of the local correction matrix with the discrete kernel renders the exact 
kernel. 
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The idea of the local corrections is to introduce a new quadrature rule that is specialized to the 
convolutional integral and that is exact for a specific class of functions. Assume that the current 
density J is expanded into a set of known basis functions f k ( r ) that are distributed over Q. Then, 


from (2.5), define the operator: 

L m J k (r n ) = J G {r m - r)f k (r')ds' -^(o n G m J k (r n ) 

»=i n "=i 

n*m 


( 2 . 6 ) 


Expressing this for K basis functions f k { r ) leads to a linear system of equations on the left-hand- 
side of (2.6), which is used to solve for the the m-th row of the local correction, L m n . The right- 

hand-side of (2.6) corresponds to the residual error of the discrete approximation of the 
convolution relative to the exact convolution of the kernel with f k ( r ). The linear system of 
equations is solved using a direct LU factorization if K = N or singular value decomposition if 
K*N . 

For problems that are sufficiently large, L m n is a sparse matrix. This arises since away from 
the singularity G(r m - r' t ) becomes a smooth function. Subsequently, in this region numerical 
quadrature approximates (2.1) to high-order. Thus, as \r m -r'\ becomes sufficiently large, the 
residual error on the right-hand-side of (2.6) tends to zero and entries of L mn become negligibly 

small. 

Once L mn is constructed, (2.3)-(2.5) are combined, leading to the linear system of equations: 

4'K.) = Io>.(G„+i„)JK) < 2 - 7 ) 


which is used to solve for the unknown current density vector. 

The Nystrom method has the distinct advantages that: 

• The discrete kernel G m , n is computed simply from one point evaluations of the kernel. 

• The local correction matrix is sparse, and only needs to be computed in the near vicinity 
of the observation point. Thus, constructing L m n requires 0(N) operations. 

• Increasing the convergence order simply requires an increase in the number of quadrature 
points. 

• The method provides a smooth approximation for the currents over the surface Q that is 
represented by an expansion of the functional basis of the underlying quadrature rule. 
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3. Application of the LCN Method to the EFIE 


3.1 The Nystrom Formulation of the EFIE 

Consider the electromagnetic scattering by a three-dimensional PEC object situated in a free 
space. The relationship between the equivalent currents on the surface of the object and the 
scattered electric field can be expressed as: 

E s = -jk oVo f ?>' - j Y V f ( V G(r,?') • J(f'))ds' (3.1) 

n ” Q 

O & 


where G(r, r') is the ffee-space Green's function: 


G(r,r') 


e -jk\r-r'\ 

47i|r -r'\ 


(3.2) 


Initially, it is assumed that the surface S is discretized using curvilinear patches that represents 
the surface curvature to high-order. It is assumed that the surface of the patch can be uniquely 
described by a two-dimensional curvilinear space (u l ,u 2 \. The unitary vectors and a 2 


defined for this space are described as [7]: 

S , = —t (3-3) 

du l 

where r = r(u l ,u 2 ). By definition, a x and a 2 are independent vectors that are tangential to the 
surface. Next, the coefficients of the metric tensor are defined as 


Si,j = a i • a j ■ 


(3.4) 


The Jacobian is the square root of determinant of the metric tensor. Hence, the Jacobian for a 
surface integration is 

■Is = V£l,l£2,2 - Si,2 • ( 3-5 ) 

The PEC body is assumed to be illuminated by an incident wave, defined by the electric field 
E mc . The resulting total tangential electric field must be zero on the PEC surface. Then, given 
fitotal _ ginc + ^ p j) the total tangential electric field is constrained to zero on the PEC 

surface, resulting in the electric field integral equation (EFIE): 

-Sj • E inc = -jk oVo J (3, • J(?'))G(f,f')ds' (v G(r,r') ■ J(f'))is' (3.6) 

s 0 s 

where a i (i = 1,2) are the unitary vectors. 
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The focus of this section is the Nystrom solution of (3.3) for the equivalent surface current 
density J(r'). Following the scheme presented in Section 2, first the underlying Nystrom 

formulation must be derived as in (2.3). To this end, a numerical quadrature scheme is 
introduced over each curvilinear patch. Assuming that there are P patches, and N Qp quadrature 

points on the p th patch. Furthermore, the vector surface current density J(r') is projected onto 

the two independent unitary vectors. Hence, it is approximated as: 

J(r,'i) = ln n V A- (3.7) 

where J, and J 0 are constant coefficients affiliated with the abscissa point , at and a' 

1» 4 n n 

are the unitary vectors evaluated at r r ', and is the Jacobian evaluated at r' n . Finally, 

applying numerical quadrature to (3.6) and using (3.7) results in the discrete equation: 

p N i P J, at +J 9 a' ) 

-a t • E mC (r m ) * -jk oVo " rr~ ^ G fmfnWnUn 


p~l n~\ 

P N a „ 


-\G{r m X)4i u 


„ P q P i (A at +J 0 a' )) _ 

-WG(44)- n n rr n - ¥n“n 

k o m 1 l U JJ 

where (r n ,cj n ) are the abscissa and weight pair for the p th patch, and for i = 1,2. 

The expression in (3.8) leads to the discrete kernel G m ,„, which can be expressed as a linear 
system of equations: 


\n ' ( r ^)] _ 


G m,n i i2 lf J l n 


E inc (7 m )\ 


m,n u 


"2,2 V 


In the second term of (3.8), there is a double gradient operation. Since the surface current 
density is expressed simply as a vector with unknown constant coefficients, the derivatives of the 
gradients must be restricted to operate on the free-space Green's function. This operation must 
then be derived to explicitly express the entries of the discrete kernel. 

In order to understand more concretely specific terms of the discrete kernel, the free-space 
Green's function is expressed in terms of its real and imaginary parts. Specifically, (3.2) is 


written as: 


p j 1 [~cos(£|r -r'l) sin(£|r —r 'I) 

G(r,r') = G R (f,n-jG'(F,n = -- 3-1 "- -j 

4rc \r-r’\ \r —r \ 


(3.10) 


It is immediately seen that G R is singular, whereas G 1 is a regular function (i.e., G 1 and all 
derivatives of G 1 are bounded at vanishing separation of r and r '). 
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Using this decomposition of the Green's function, the V(V G(r,r) ■ J(r )j term can be 
evaluated separately for G R and G 1 . It can be shown that: 


i7 -v((vGV,r'))••/(?')) = -a r J(r')k 


3 sin kR cos kR 
(kR) 2 (kR) 3 


_ _ _ - „cos kR „sin kR cos kR 

+a t -(r -r )J(r )-(r -r )k 3- - + 3 - — - -r 

[ (kR) 5 (kR) 4 (&?) J 

« j _\ - _ \ _ - ,3 cos kR sin kR 

v A ; L ( kR ) ( kR ) w J 

where R =\r - r'\ . In deriving (3.1 la) and (3.1 lb), the vector identity is used where: 

. J(?-?■) a, ■}(?') (*r(f-np(r'Hr-n) 

‘ V R ) R R 3 


(3.11a) 


(3.11b) 


(3.12) 


The expressions from (3.10) and (3.11a,b) can then be used within (3.8). To this end, the 
discrete kernel in (3.9) is expressed as: 

G m ,n. . = Gm,n . + &m,n ., • (3.13) 

’ h3 ’ h3 h3 

where, G R n arises from G R , and G l m ^, from G r . Subsequently, 

.7] 0 ,2/-> -t \ coskR m,n s ^ nkR m,n cos 

i kRm j 

[l > V *™> n I (3. 14 a) 

, 74/-* B \/-/ /-> ->/ o CQS kR vi,n Q s * n kR rn;n cos kR m,n 

\im ' Rm t n ) \j n ^™ rn G (bR \ 5 (bR l 4 (kR. _ ) 3 


sin kRynn 

cos kRjnn 

(«W) 4 

(kRmn) 
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= ^o 
4t r n 


k 


K. - 


Jn 


sin kR^ 

sin kR 

kRjn.n 


cos kR^fi 


[kRm,n ) 


+ 


• Rm,n){^j n ' Rm,nj 


sinfci^,, 
sin kRfn^ ^w,n 


cos kRfftj 


kRm,n 


(kR m ,n ) 


{kRm,n ) 


(3.14b) 


where Rm^ — ( r m r n ) snd 

The expression in (3.14a) is clearly hypersingular in the limit that R^ n —* 0. Subsequently, 

this term must not only be locally corrected for a well-posed Nystrom formulation - this term 
must also be specially treated such that the numerical integration required for the exact kernel is 
convergent. On the other hand, as expected, the integrand of (3.14b) is regular and hence is 
integrable to high-precision using a numerical quadrature. This is observed by performing a 
Taylor series expansion of the trigonometric terms and observing their value in the limit that R 
tends to zero: 
sin kR' 


lim 


lim 
R —>0 


kR 


sin kR 
kR 


= 1 


cos kR 


(3.15a) 


(kRY 


kR- 


(kRf 
3! 


+ 


kR 


( 1 - 


{kRf 

2! 


+ •••) 


(kRy 


1 

3 


(3.15b) 
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lim 

0 


sin kR 

sin kR 3 kR 
kR 


— cos kR 


(kR) 


(kR) 1 


kR JkKf_AkRf 

kK 3! + 5] 


(3.15c) 


(kR) 2 (kR) 4 

[ 2! 4! 



15 


Subsequently, from (3.14b) and (3.15a) - (3.15c), it is concluded that: 


lim < n . . = ^ J, u n 

hi 47T ^ 


m-+n 


[a- • a' 

\ hn Jn 


-1 + ' 


• Rm,n){^j n 4 Rm,n ) 15 


= 9slJ' u 
4t r A. n 


1,2/-* "+/ 

- k a • a • 

3 \ Ji 


(3.16) 


^771 


In summary, (3.8) - (3.16) present the underlying Nystrom formulation for the EFIE. It was 
observed that since G 1 is regular, then the contribution due to G l in the discrete integrand of (3.8) 
is regular and can be evaluated to high-order via numerical quadrature. Subsequently, these 
terms do not require local corrections. On the other hand, G R is singular, and the terms of 
Gm n actually exhibit hypersingularities. Subsequently, these terms must be locally corrected 

’ hj 

to realize both a well-posed formulation and a formulation that yields high-order convergence. 
Furthermore, due to the hypersingularity of the integrand, this term will have to be manipulated 
such that the numerical integration is numerically tractable. 
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3.2 Locally Correcting G R 

3.2.1 Patch Based Local Corrections 

The local corrections of the discrete kernel will be performed in a manner similar to 

5 

that outlined by (2.6). Specifically, specialized quadrature rules will be developed such that the 
evaluation of (3.8) is exact for a given set of basis functions. Interestingly, only the contribution 
due to G r will be locally corrected. 

As outlined in Section 2, the local corrections are performed by first introducing a set of 
suitable basis functions that are distributed over S. This basis set will actually serve as a sample 
set of currents distributed on S. The field radiated by these currents will be calculated at the field 
point r m on S. This will be calculated in two ways. The first is using the discrete kernel in (3.8) 

(with G„ im =0). The second will be using an exact integration of (3.6) for the kernel G R . 

’ h3 

Since this expression cannot be evaluated in closed form, the exact integration is performed 
using an adaptive quadrature to a specified precision. The residual representing the difference 
between the exact integration and the discrete kernel is then computed for each basis function. 
This residual term, which appears on the right-hand-side of (2.6), will be used as the forcing 
function for the local correction terms. 

Two classes of local correction operators can be used to correct the discrete kernel. The first 
is what is referred to as patch-based correction. The second is referred to as entire-domain 
correction. Patch based local corrections imply that the new quadrature rule represented by the 
local corrections is limited to the support of a patch. Thus, for each field point, a new quadrature 
rule is established for each patch. For the entire-domain approach, a new quadrature rule is 
established for a sufficiently large area surrounding the field point. 

Patch-based local corrections restricts the support of the basis functions to that of a patch. 
Subsequently, this limits the support of the quadrature rule derived for each local correction. 
The advantage of the patch-based local correction is that the linear systems of equations needed 
to solve for the new quadrature weights are small and generally well conditioned. It also allows 
the use of orthogonal functions that lead to orthogonal matrices that are analytically solved. 
Another advantage is that the corrections rely on a local curvilinear coordinate system rather than 
a global coordinate system. This provides for a more robust solution for three-dimensional 
surfaces. The disadvantage is that this scheme can lead to asymmetries in the final impedance 
matrix, and if one is not careful can result in current discontinuities in the final solution. 

On the other hand, entire-domain based local corrections utilize basis functions that are 
distributed over the entire surface S. The support of integration is then limited to being 
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sufficiently large to reduce the residual error below a tolerable level of accuracy. One advantage 
of such a scheme is that current continuity is enforced by the basis functions. However, for 
three-dimensional surfaces, such a scheme is often limited to canonical geometries for which a 
global curvilinear coordinate system is easily derived. Furthermore, entire-domain schemes lead 
to large linear systems that are often very poorly conditioned. 

In this report, the focus is limited to patch-based correction schemes. Entire domain schemes 
can be derived in similar manners - but without the bounded support of the integrations. The 
details of the corrections are provided below. 

Each local correction is associated with a single field point r m . Thus, for a patch based 
scheme, only the row entries of G m , n associated with a single patch are corrected at a time. A 
two-dimensional basis set is introduced over each patch. Within the curvilinear coordinate 
system, the vector basis is introduced as: 



where k represents the order of the testing function, and the a' are the unitary vectors evaluated 


at 



It is noted that the basis is normalized by the Jacobian such that the current basis is 


divergenceless. The vector current basis will serve as a sample current distributed over the 
patch. Subsequently, the field radiated by the current is calculated at r m using exact integration 

and the discrete kernel. The local correction establishes a new quadrature rule which when 
operating on the current basis yields the residual error between the exact and discrete 
integrations. This is actually applied independently to each quadrant of the discrete kernel in 
(3.9). For example, consider quadrant (/, j). The quadrature rule operating on the &-th basis 
function is expressed as: 

= -XoVof (*„' S')G S (4,?') 

n=l s 





i 


vg k (4,?')' 



/ 


(3.18) 


- E <», ,//KK 

n =1 
n^m 


where, r m is the field point, and r' n and u n are the abscissa and weights of the underlying 
quadrature rule and the patch being corrected, and 7 n are the new quadrature weights to be 
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computed. The first two terms on the right-hand-side of (3.18) correspond to the exact 
integration of the convolution of the kernel with the basis functions /* (r), and the last term the 

numerical quadrature approximation of the convolution with defined in (3.14a). This 

procedure is repeated for K basis functions, leading to a linear system of equations: 



(3.19) 


where R ^ is the residual error resulting from the right-hand-side of (3.18) for the fc h basis 
function /*(r). This is then solved for the vector -f n . Subsequently, the product u; n 7 n is then 


inserted in the appropriate entries of the m-th row of the local correction matrix in the (i, j ) 
quadrant. 

Once the local corrections are performed for the m-th row, then the impedance matrix is 
updated. Specifically, from (2.5) each quadrant is corrected as: 


G. 


m,n M 


= <7 + L 


7 / 1 , 71 . . 


(3.20) 


3.2.2 Evaluating the Self Patch 


In the preceding section, the residual term on the right-hand-side of (3.19) requires the 
evaluation of the exact integral operator operating on the basis function. When the field point 
lies on the same patch as the basis function, the integrand will be hypersingular, as described in 
Section 2. Subsequently, the integral operator has to be modified such that the integral is 
numerically tractable. Specifically, the hypersingular terms arise from the gradient divergence of 
the real part of the Green's function. It is this term that will be manipulated. To this end, this 
term is first rewritten as: 


J dj • v((vG*(r,r'))- •/(?'))*' = - J a t • v((v'GV,F''))• Jir'^ds' 

S S 

= -jj(r’)-Vl l {d r VG R (r,P')y s ' 

S 


(3.21) 


where WG R (r,r') = -V'G R (r,r'), and V| is tangential to the surface. Next, utilizing the 
vector identity V • = (j>V • A + A • V<j>, the right-hand-side of (3.21) is rewritten as 
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(3.22) 


It is then realized that the first term on the right-hand-side of (3.22) can be rewritten using the 
divergence theorem for open surfaces [2] as: 



a t -VG R (r,r ')j ds' = <£(e'-/(r'))| 


ai • VG 


R 



(3.23) 


where C is the closed contour bounding the open surface S, and e' is the outward normal to the 
contour that is also tangential to the surface (i.e., e'dl' = dl' x«). The second-term on the right- 
hand-side of (3.22) is rewritten as: 

J V| J(F')(«/ •VG R (r,r')jds' = j* VG^(r,r')-[a,-V| • J(f')]<is' (3.24) 

5 5 


This is then rewritten as: 

= J VG R (r,r’) -[a/Vfi • J(r ')- kfs' + J VG R (r,r')■ K'ds' (3.25) 

5 5 


where, K' is a vector defined by 


K' = al 


Vfi -J{r') 


r =r 


ai 


Vi 7 '" 


The constant % is chosen such that at the singular point K' = a t V[| • J(r'). Thus, 



(3.26) 


(3.27) 


Thus, choosing K' according to (3.26) and (3.27), the integrand of the first integral in (3.25) 
becomes regular. This occurs because the cancellation of a,-V| • J(r') cancels out one of two 

poles of VG R (r,r'). This leaves an integrand with a single pole that can be integrated 


numerically. 

The second term in (3.25) must still be manipulated into an integrable form. To this end: 


J VG R (r,r')- K'ds' = -J V'G R (r,r')-^Lxds' 


s 

— 


— vigV.?') 


Vi 7 


ds' 


(3.28) 


= “ X 


l 


Vfi 


a>i G R (r,r’) 


VF 




:r\vn. . a i 


\ds f 
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It can be shown that V| 


a i — = 0. Thus, applying the open surface divergence theorem: 


f VG* (r,r')- K'ds' = -xfe' • -% G* (P.r')dl' 

i c & 


(3.29) 


Finally, from (3.21)-(3.29), 

o. • Vj (vG R (r,r') • J(r'))ds' = +f VG R (r,f') ■ ^V' • J(r') - K 


ds' 


- §(e' ■J(?'))(a i -VG ft (r,r'))dl’ 
G 

n <9 


(3.30) 


The expression in (3.30) can now be used to evaluate the right-hand-side of (3.18). To this 


end, define: 
FtR 


hJ 


-jk 0 Vo f( S i ■ a'j)G R (r m ,r')fj (r')ds' 
s 

- j^a t ?')■«'// (?'))*' 

c V 


(3.31) 


where the second term on the right-hand-side of (3.31) is to be replaced with (3.30). 

The above integrals are performed using an adaptive quadrature. The leading singularity is 
1/R, which is numerically tractable. To accelerate the integration, a Duffy transform is used [8]. 
To this end, the patch is first subdivided into triangles with sharing a common vertex at r m . The 

Duffy transform simply defines the triangle as a quadrilateral with two vertices shared at the 
common vertex at r m . Subsequently, the Jacobian goes to zero at this vertex. This greatly 

accelerates the adaptive quadrature. 
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4. Application of the LCN Method to the MFIE 
4.1 The MFIE 

An integro-differential operator magnetic field integral equation (MFIE) can be used to solve 
for the induced equivalent surface currents providing the PEC body is a closed surface. The 
MFIE is applied on the surface of the body through the simple boundary relationship: 

where a n is the outward normal, S + indicates the closed surface just outside of S, and S~ the 
closed surface just inside of S. Next, a test vector a % that is tangential to S is introduced. 

Subsequently, (4.1) is modified as: 


®j ’ Q'n ^ ^ H 




—a ■ H tot , = a,- dr, xj 


+ ~ “i l“» 


Finally, letting H tot = H s + H mc , and realizing that 
H s (r)= -f VxG(r,r')j(r')ds' 


leads to the MFIE: 


-Sj • H inc (?) = Sj ■ (a„ x J (?)) + a, • f V x G(r,f’)J (r ')*' 


where, r € S + . 

The surface current density is again assumed to be expanded in terms of two independent 
vectors chosen to be the unitary vectors defining the surface as: 

jr(f) = (4.5) 

where Jg is the Jacobian and a, t (i = 1,2) are the unitary vectors at r . It is noted that from 
(3.4) and (3.5), it can be shown that 

Jg=\a 1 xd 2 \. (4.6) 

Subsequently, from (4.5) and (4.6), the first term on the right-hand-side of (4.4) is evaluated 


a. -(d n xJ (r )) = n t 


-J v i = 1 
J -,, i — 2 
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Next, the integral on the right-hand-side of (4.4) is analyzed. To this end, it is realized that. 
Vx(j (r') G (f, f ')) = G (r , r')V x J (r') - J (f ') V G (r, r') 

= -j(r')VG(r,r') . (4.8) 

= j(r')V'G(r,r') 


where V' implies the gradient in the primed coordinate system. Next, it is seen that: 

—► 

V' G (?, ?') = -1 ^ \g r (R) + jG’ (J8)]. 


(4.9) 


where R = r — r',R = 


, and G R and G 1 are defined in (3.10). This is explicitly written as: 


V'G (?,?') = ~ R\K R (R) + ]K\R) . 

V ’ 47T L j 


(4.10) 


where, 


%. . kR sin ( kR ) + cos ( kR) 


K (R) = 


(kRf 


(4.11) 


and 


K’(R) = 


cos 


(kR) 


sm 


(kR) 


kR 


(kR) z 


(4.12) 


Ag ain , it is seen that the contribution from the imaginary part of the Green's function is regular, 
and that from the real part is singular. 

The expressions from (4.10) through (4.12) can then be used within (4.4). Upon inspection, it 
becomes obvious that the integral over the real term K R (R) requires the principal value integral. 
Specifically, it is assumed that r is on S + . However, K r (R) is hypersingular on S + . Thus, 
this integral must be manipulated. To this end, define r is on S + as: 

r = rg + a n e . (4.13) 

where r g is on S, a n is the unit normal to S at r s , and e is the projection between S + and S 


(e >0). Then, one can also define: 

R = r —r' = r s —r' + a n e = R s + a n e . 


(4.14) 


where in the limit that e —> 0, 


R 



Then, from (4.7), (4.10), and (4.14), (4.4) can be 


rewritten as: 
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-h-« inC (?) = "i + 


47T 


J R s ■ (s i xJ(r’))(K R (R)-jK' (iS))V7*' 

S 


v k 

+ lim — 
£—►() 47T 


f d n e-(a t x J^K^-jK^R^ds' 

s e 


(4.15) 


where k is given by (4.7), and the second term on the right-hand-side is a principal value 
integral over the small surface S £ encircling the field point. Focusing on the first term on the 
right-hand-side of (4.15), the integrand involving K ! is regular. The integrand involving K R is 
singular and integrable. Specifically, from (4.11), K R has a 1 / R 3 singularity. However, in the 
limit that R -> 0, both R s -> 0 and the dot product R s -a t xJ (r') -> 0. The latter condition 

-4 

arises since the cross product is normal to the surface at the field point, and R s is tangential. 
Thus, the integrand will have a 1 / R singularity, which is integrable. 

Next, the principal value integral is treated. Initially, since K' is regular and bounded as 
R -+ 0, the integral of this term will tend to 0 as e -» 0. However, based on the previous 
arguments, since the dot product a n ■ (a- x J (r')j is not zero in the limit R —* 0, the integrand 


is hypersingular. Let us consider this term further. For simplification, consider the case when 
a % = a 1 . Given J (f') defined by (4.5), then, it is seen that: 


lim cL 




' 2 - 


(4.16) 


Since the contribution due to K 1 is zero, the principal value integral in (4.15) can be rewritten 
using (4.11) and (4.16) as: 


PVI= lim— eJ 0 -f 
£-►0 47r 

Se 


kRsin(kR) + cos (kR) 


(kRf 


Jg'ds' 


(4.17) 


In the region near the field point, R is assumed to be very small. Subsequently, small argument 
aproximations for the trigonometric functions can be applied, where: 


sm(kR)^kR + 0(ikR) 3] j 
cos (kR) ~ 1 T o[{kR)^ j 


(4.18) 


Next, the surface S £ is approximated by a small disk of radius A p centered about the field 
point. Subsequently, the principal value integral is written as: 
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b 6 

PVI= lim— eJ 0 

£__>0 47r ^ 


77 ^ 

p==0 (p—O 


+ 


(kRy (kR) ) 


pd(j)dp 


Next, let R = \]p 2 + £ 2 , which leads to: 


k 6 

PVI = lim — eJ 0 

e —>0 2 2 


k" 


A p 

I 

P =o 




W p 2 +£ 2 fc 3 (p 2 + e 2 ^ 


pdp 


= lim — eJ. 


e —>0 2 
3 


£ 2 -£yj 

^A p 2 +£ 2 ,-£ + •> 

jAp 2 + £ 2 

k k 3 f 

A p 2 +£ 2 


k 


J n 
2 2 


0 + 


Jo 


(4.19) 


(4.20) 


Similarly, for i = 2, it can be shown that the PVI = -- J 1 . Finally, from (4.20) and (4.15), the 

2 


MFIE is expressed as: 


«.• r 


-a, ■H’ nc (r) = ^ + yfM s -(a i xJ (?'))(k R (R) - jK 1 (/?)) ,/?*' 

^ —t/l 

5 


(4.21) 


4.2 The Locally Corrected Nystrom Formulation of the MFIE 

The MFIE will also be analyzed using the locally corrected Nystrom scheme. To this end, the 
surface S is discretized into P patches. As for the EFIE, it is assumed that the patches are 
curvilinear and provide a high-order description of the surface geometry. Subsequently, (4.21) is 
approximated by numerical quadrature by introducing a quadrature rule for each patch. Then, 
sampling the field point at abscissa point r m , the Nystrom approximation of (4.21) is: 


K- 


~\i ' H mC = ~2 ^ m ' n 

3 p N q ^ (4.22) 

+ — ^ ^ Rm y n * x J (^n )) [f R (^m,n ) — 3% (^m,n )) u n yfdn 
4 %=ln=l 


where, ^ =r m -r n , i^ n = R m>n , and 8 m ^ n is the Kronecker delta function. 
Thus, a discrete kernel is defined as 


p _ pR i pi 

■ — - ' ^m.n. . 

’ ’ h3 i,3 


(4.23) 


where, the m,n ih entries in quadrant (i,j) of the discrete kernel is: 
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G k\j (% n x n 'j n )( KR { R m,n )) 

G k\j = JZ Rm,n ‘ (\ n X X ) ( Rl (^n,» )) 


(4.24a) 


(4.24b) 


where, o is a unitary vector evaluated at the field point r m and a'- is a unitary vector 

’ hn Jn 

evaluated at the source point r n . 

In the limit that m = n (i.e., -* 0) K 1 is bounded, as seen from (4.12). Furthermore, in 

this limit, the dot product is zero, and R mn Subsequently, G^ro- =0. On the other 

hand, G5, m is infinite in value due to the singularity of K R . Thus, G* n . . must be loca,1 y 

’ Ifl,III 


corrected. 

The local corrections can be carried out in the same fashion as outlined in Section 3.2.1. One 
advantage of the MFIE over the EFIE is that the level of singularity of G^ m is simply 1/R. 

Subsequently, the exact integration in (4.21) can be carried out numerically in an efficient 
manner using adaptive quadrature and the Duffy transform [8]. 
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5. Numerical Validation of the Nystrom Solution 

The locally corrected Nystrom (LCN) methods presented in this paper report been applied 
successfully to the electromagnetic scattering analysis of PEC bodies. This section summarizes 
some of these results. The focus of this study is on the accuracy and convergence characteristics 
of LCN method when applied to the computation of the RADAR cross section (RCS) of PEC 
objects. Initially smooth bodies are studied. Then, the method is applied to open structures with 
edge singularities. 

The error convergence of the LCN method will specifically be presented by studying the root 
mean square error (RMS) in the RCS as the level of discretization is modified. The RMS error is 
computed relative to a known reference solution which is either obtained analytically for 
canonical geometries or numerically through a LCN solution with a very high level of 
discretization. Specifically, given the RCS computed over a sphere in the far field, the RMS 
error is estimated as: 


RMS Error = 


27r 7T 

J J | RCS ref (9,<f>) - RCS lcn {9, <f>) | 2 sin 9d9d(p 

4>^oe=o __ 

27T 7T 

J J \RCS rcf (9,4>) | 2 sin 9d9d(f> 

( p —0 0=0 


(5.1) 


where 9 and (j) are the observation angles, RCS r ^(9,(j)) is the reference RCS computed 
analytically (if possible) or numerically to high precision, and RCS LCN (9,(j >) is the RCS 


computed via the LCN method. 

The first example illustrated is that of a PEC sphere. The study of this problem is of interest 
because of its smooth geometry and the fact that the RCS can be calculated analytically using a 
Mie series solution. The sphere was discretized using curvilinear quadrilateral patches that 
exactly modeled the curvature of the sphere. At the poles of the sphere, the quadrilateral 
essentially becomes triangular in shape as two vertices will coincide. The quadrature rule for 
each sphere was then determined via a product rule of one-dimensional Gauss-Legendre 
quadratures. Orthogonal patch based testing functions were used for the local corrections. The 
number of testing functions used was equal to the number of quadrature points, resulting in an 
orthogonal test matrix. 
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Fig. 5.1 RCS of a PEC sphere versus $ computed in the 4> = 90° plane when illuminated by an x- 
directed electric field propagating along the negative z-direction with radius defined by k G a = 1. 
Results from the locally corrected Nystrom method with a total number of unknowns equal to 64 
and 400 and an exact MIE series solution are illustrated. 
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Fig. 5.2 RMS error in the RCS for the LCN solution of the MFIE for the scattering by the k Q a = 1 PEC 
sphere as the quadrature order is increase, (a) RMS Error relative to the Mie series solution, (b) 
CPUs required to complete the computation on a MIPS R10,000 processor. 
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The sphere is assumed to have a radius a such that k 0 a = 1, where k 0 is the free space 
wavenumber. The sphere was illuminated by a plane wave traveling in the negative z-direction. 
The sphere was discretized using 8 quadrilateral patches (i.e., one patch per octant). The RCS 
was then computed as the quadrature order was increased. Figure 5.1 illustrates the RCS 
computed by the LCN method for two levels of discretization. The first used 4 quadrature points 
per patch (2 along each linear direction), which resulted in a total of 64 unknowns. The second 
assumed 25 quadrature points per patch (5 along each linear direction), resulting in a total of 400 
unkowns. These results are compared to the exact Mie series solution. 

The case with 4 quadrature points per patch has an error of roughly 0.2 dB in the backscatter 
region (0 = 0). The 25-quadrature point per patch case is indiscernible from the exact solution 
on this graph. 

The RMS error in the RCS was then calculated as a function of the number of degrees of 
freedom. These results are illustrated in Fig. 5.2(a), while the CPU time recorded for these 
simulations on a MIPS RIO,000 processor is illustrated in Fig. 5.2(b). For these calculations, the 
degrees of freedom were simply increased by increasing the quadrature order per patch. 
Subsequently, on a log-log scale, the RMS error recorded in Fig. 5.2(a) is increasing in slope as 
the quadrature order is increased. This clearly shows that the Nystrom scheme realizes 
exponential convergence. The case of N = 256 realizes 4 digits of accuracy, requiring only 8 
CPU-seconds to compute, which is more than sufficient for engineering accuracy. 

As demonstrated by the previous example, the LCN method yields high-order convergence for 
the RCS computation of closed smooth bodies (More extensive studies of this can be found in 
[2]). The question that arises is how applicable is the LCN for analyzing the RCS of surfaces 
with edge singularities. Specifically, when modeling geometries with edge singularities, a high- 
order solution can only be achieved if the surface current density can be represented to high 
order by the underlying quadrature rule. Thus, standard Gauss-Legendre quadrature would not 
provide high-order convergence, unless the discretization were refined in an appropriate manner 
near the edge. This would lead to an unnecessary increase in the number of degrees of freedom 
and the condition number of the impedance matrix. A more viable alternative is to introduce 
specialized quadrature rules that implicitly integrate the singular behavior to high order. Such 
quadrature rules have been successfully implemented for sharp edges on thin sheets in two and 
three dimensions [3, 9, 10]. It was shown for two-dimensional geometries, that a high-order 
representation can be realized using specialized quadrature rules [3]. As an example, consider 
the edge currents of a flat plate. The normal edge current are oc x l/2 , and the tangential currents 
are oc x ~ l/2 , where jc is the distance from the edge. In three-dimensions this current can be 
integrated to high order using a Gauss-Jacobi quadrature rule [11]. 
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The local corrections require the choice of a suitable basis f k (u\u 2 ) to be distributed over Q. 

In this work, orthogonal subdomain functions representing the underlying quadrature rules are 
assumed. Specifically, for a quadrilateral cell, f k (u l ,u 2 ) is equal to the product of two one¬ 


dimensional Jacobi polynomials: 

t(u'y)=p^ (u')p^ (u>) 


(17) 


where P^(u) is the k-th order Jacobi polynomial, (u\u 2 ) are the parametric coordinates for the 

cell with each defined over the range (-1,+1). For smooth patches without an edge singularity, a 
and 0 are chosen to be 0, realizing that (u) is equal to the k-th order Legendre Polynomial. 

This choice of functions has a number of advantages. Initially, the local corrections can be 
efficiently performed on a cell by cell basis. Secondly, if appropriately scaled, the matrix 
resulting on the left-hand-side of (6) is orthogonal and the matrix can be diagonalized 
analytically. Finally, with the use of Jacobi polynomials, a and/5 can be appropriately selected 

such that singular currents are represented to high-order. 

As an example, consider the finite strip array illuminated by a uniform plane wave illustrated 
in Fig. 5.3. The bistatic RCS for this array was computed using both the LCN and a 

commercially available method of moment (MoM) simulator (Zeland IE3D). For the LCN 
simulation, each strip was discretized into 5 rectangular cells, each being 0.25 A 0 x 2A 0 in 


dimension. The discretization was then refined by increasing the quadrature order. The MoM 
solution employed rectangular cells with roof top basis functions. Special edge cells were also 
included to improve the edge-current representation. The MoM discretization was refined by 
increasing the number of patches. The strip array was illuminated by a uniform plane wave 
incident from the broad side as illustrated in Fig. 5.3. The RCS in the 0 = 0° plane as computed 

by both a MoM simulation with over 3000 unknowns and the LCN method withl080 unknowns 
is illustrated in Fig. 5.4. Excellent agreement is realized in this cross-section. 

The root-mean-square (RMS) error in the RCS was then calculated for each technique as a 
function of the discretization level. These results are presented in Fig. 5.4. The error 
convergence realized by the MoM solution is only first-order. In fact, to realize 1 % error for 
this problem, the LCN method requires roughly one-fourth of the computational resources as the 
low-order MoM solution. Obviously for errors < 1 %, the LCN requires only a modest increase 
in resources. Whereas, it is unreasonable to expect better than 1 % accuracy from the low-order 
MoM solution. This demonstrates the specific advantage of the high-order scheme that increased 
orders of accuracy can be realized with only modest increase in the computational resources. 
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Fig. 5.3 Finite array of thin PEC strips in free space Fig. 5.4 Bistatic RCS of the finite strip 

array computed via the Nystrom 
method. 



Fig. 5.5 RMS Error in the RCS versus N for 
both MoM and Nystrom solutions. 


Fig. 5.6 CPU time (in seconds) versus N for the MoM 
and the Nystrom solutions. 
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The CPU times for the simulations in Fig. 5.5 were also recorded. For this calculation, the 
LCN solution was performed on a MIPS R10000 180 MHz processor with 1 Mb cache. The 
MoM solution was performed on a 350 MHz Pentium III PC. Thus, one cannot compare directly 
the CPU times. However, one can compare the rate of increase of the CPU time with the rate of 
error convergence. Specifically, the MoM simulation is linearly converging as the CPU time is 
increasing as 0(N 3 ). On the other hand, the LCN simulation is exponentially converging as the 

CPU time is increasing as 0(iV 2 ) as the matrix fill dominates the computation as compared to 
the factorization. 

The final example illustrated in this report is the scattering by a PEC hemisphere. The 
hemisphere has a radius defined again as k B a = 1. This problem poses a particular challenge 

that the geometry is curved and it contains an edge singularity. The structure is discretized using 
curvilinear quadrilateral patches as used for the spherical scatterer. As with the strip study, 
Gauss-Jacobi quadrature is defined for patches containing an edge singularity and Jacobi 
polynomials are used for the basis functions for these patches as well. It is realized that since 
the surface is non-planar, that the singular behavior is not modeled exactly by the quadrature 
rule. Neverthless, it is anticipated that the rate of convergence will be still superior to that 
obtained using straight Gauss-Legendre quadrature. 

A simplified geometry for the hemisphere is illustrated in Fig. 5.7. It is noted that the 
discretization and shading are an artifact of the imaging software and do not represent any 
physical meaning. The hemisphere was illuminated by an incident plane wave as illustrated in 
Fig. 5.7. The RCS computed in the 0 = 0° plane is illustrated in Fig. 5.8. The RCS was 

computed both using a MoM solution (Zeland IE3D) and the LCN method. The convergence in 
the RCS as a function of the discretization is illustrated in Fig. 5.9. For this study, the base 
discretization consisted of eight curvilinear patches that exactly modeled the surface curvature. 
To this end, there were two patches along the theta-direction (or zenith), and four along the phi- 
direction (or azimuthal). The number of quadrature points was then increased from 4 points per 
patch to 49 points per patch. The reference solution for the RMS error calculation was a more 
finely discretized solution with 12 patches and 81 points per patch. 

Observing Fig. 5.9, it is seen that the RMS error is not demonstrating high-order convergence. 
The convergence is still more rapid than that expected for a low order method. However, it is 
not exponentially converging as the previous cases. It is anticipated that this is due to the fact 
that the edge singularity is not being modeled exactly by the Jacobi polynomial. However, since 
the edge behavior is better approximated than using a Gauss-Legendre quadrature rule, the 
convergence is still improved. 
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Fig. 5.7 Plane wave illumination of a PEC hemisphere 

ka = 1 . 


Fig. 5.8 RCS of the PEC Hemisphere versus Om the 
cj> = 0 plane. 
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Fig. 5.9 Error convergence of the LCN solution of the 
scattering by the PEC hemisphere. 


Fig. 5.10 CPU time for the calculations in Fig. 5.9. 


DAA G55-98-1-0460 


27 


The University of Kentucky 

























STIR: High-Order Sparse Matrix Nystrom Algorithrn for the Analysis of Planar Microwave Circuits - F/na/ Report 


6. Multilayered Green's function Computation 

6.1 The MPIE 

The Nystrom scheme presented in the previous sections has been developed for the tree-space 
Green's function. This is a necessary step in the development of this program. However, the 
ultimate goal is to apply the Nystrom method to the solution of printed circuits in a layered 
media. To this end, the Nystrom scheme will be applied to the mixed potential integral equation 
(MPIE) form of the EFIE. Specifically, (3.6) is rewritten as: 

-a { • E inc = —jk 0 %Ja,-5 A (?,r')-J(r >'+ iy v/G,(?,/)(v J(r'))*' (6.1) 

s s 

where G. is the dyadic Green's function representing the magnetic vector potential and G q is 

A 

the scalar potential for the layered medium. 

In general, the Green's function in a layered media is not known in closed form, but is actually 
computed numerically [12]. The most computationally intensive component of this is the 
computation of Sommerfeld integrals that arise for both the magnetic vector and electric scalar 
potentials. For the Nystrom simulation, it is imperative that this integrals can be computed 
efficiently and to controllable accuracy. Traditionally, these integrals are tabulated a priori and 
when performing the integrations for the computation of the impedance matrix are computed 
using interpolation. A similar scheme is proposed here. However, since the Nystrom method is 
a higher-order scheme, the Sommerfeld integrals must be evaluated to controllable accuracy. 
Secondly, the precomputation of the Sommerfeld integrals at discrete points must also be 
performed efficiently and to controllable precision. This is addressed in the following sections. 
Initially, an efficient means of computing the Sommerfeld integrals is presented in Section 6.2. 
To this end, an appropriate contour path is introduced that avoids the highly oscillatory behavior 
of the integrand near branch point singularities and poles. In Section 6.3, an efficient and 
accurate interpolation method based on a rational function interpolation is introduced. The 
advantage of this method is that the order of interpolation can be modified to control the 
accuracy of the interpolation. Furthermore, the method reduces the number of overall points that 
need to be computed. 
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6.2 Efficient Computation of the Sommerfeld Integral 

The generalized Sommerfeld integral is given by the expression 

(6 - 2 > 

0 

where/is a function dependent on the characteristics of the layered media which is assumed to 
be singular in the region of integration and 



As an example, consider a single-layer planar dielectric substrate with relative permittivity e r 
and thickness h backed by a ground plane. The horizontal projection of the magnetic vector 
potential G at the dielectric-air interface (z = 0) and due to a source at this interface is 


expressed as: 

G A (p) = ] J o{ k pp)-^- k p dk P (6 - 5) 

0 U TE 

while the electric scalar potential G q is expressed as: 

o.M-fafcphns-*,*, < 6 - 6 > 

0 U TE U TM 

where 

N = u l +u 2 tanh (u 2 h) (6-7) 

D te -u l +u 2 coth(w 2 /?) (6.8) 

D m =e r w, +u 2 tanh(w 2 /i) (6.9) 

“,= A.< 6 - 10 > 
“ 2 = A, =(K ~ k l) in < 61 ') 


where k 2 is the wave number of the dielectric substrate, and is the free-space wave number. 
The integrals in (6.4) and (6.5) are performed on the real axis of the complex k p plane. 
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However, the path must be either deformed about the branch point at k x and the poles between 
k x and k 2 on the real axis, or a principal value need be evaluated at these points. It can be 
readily shown that D TE will not have any zeros unless 

/ >-( 6 . 12 ) 

4/i^/e r -1 

which corresponds to a relatively thick substrate. On the other hand, D TM will always have at 
least one pole on the interval of [k v k 2 ] . As the frequency is increased, additional poles will be 

introduced on this interval. 

When evaluating (6.5) or (6.6), a principal value integral must be introduced to evaluate the 
integral at any singularitie lying on the real axis. To this end, the zeroes of D m and D TE must 

be found to determine the poles of (6.5) and (6.6). Subsequently, the residues at each pole are 
then computed. Then, the principle value integral can be. It is noted that for general layered 
media, the poles must be determined numerically. However, the accuracy of the Sommerfeld 

integration depends upon the ability to accurately compute the zeroes and residues. Thus, to 
evaluate the Sommerfeld to high precision, the zeros of D m and D TE need also be determined 

to high precision. This is often computationally intensive. Furthermore, near the poles, the 
integrand becomes highly oscillatory. As a result, computing the Sommerfeld integral to high- 
precision using an adaptive quadrature scheme is also extremely expensive. 

An alternate approach to computing the Sommerfeld integration is to deform the path of 
integration off the real axis. While an infinite number of paths can be chosen, an effective 
contour consisting of a semi-elliptic path has been suggested by Gay-Balmaz and Mosig [13]. 
The semi-elliptic path suggested extends from the origin to a point on the real axis well beyond 
the zeroes, as shown in Fig. 6.1. The center of the ellipse is chosen at 

k =Mll(6.13) 
c 2 

so that the path of integration returns to the real axis well beyond any pole. 

By deforming the path in this manner, the singularities are avoided and thus the accuracy of 
the Sommerfeld integration is no longer bound to the accuracy of the prediction of the poles and 
residues. Secondly, by properly determining the minor axis of the ellipse, the highly oscillatory 
nature of the integrand is avoided. This greatly improves the efficiency of the computation. It is 
noted, though, that it is not required to use an elliptic path. In fact, a path through the saddle 
point would be ideal. However, the elliptic path has the advantage of simplicity and flexibility. 
Further, the semi-elliptic path is easily parameterized and an adaptive quadrature routine can be 
used to perform the numerical integration. 
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Fig. 6.2 Elliptic path for integrating around poles. 


Finally, the integrals in (6.5) and (6.6) are being performed by integrating along the semi¬ 
elliptical path along the interval [0,2 k c \. This is done using a simple adaptive quadrature 

scheme. Then, the method of averages [12] is used to perform the integral along the real axis 
over the interval [2A; (: ,°o]. This renders an efficient evaluation of the Sommerfeld integrals that 

can be performed to controllable precision. 

As an example, consider a single layer dielectric substrate with relative permittivity of 4.7 
backed by a ground plane. The magnetic vector and scalar potentials were computed using (6.5) 
and (6.6). In order to determine the relative accuracy of the routines, the results were compared 
to the small argument approximation of the vector and scalar potentials . Initially, the relative 
percentage of error in the magnitude of the magnetic vector potential as compared to a small 
argument approximation [14] is illustrated in Fig. 6.2 for substrate thickness of h = 10" 3 and 
10 _1 . For these computations, the adaptive quadrature requested a relative error of 10"* and 
used 25 intervals for the method of averages. It is observed that for the thinner substrate the 
small-argument approximation provides 10" 5 % accuracy (7 digits) for distances of 10' G A,or 
less. For the thicker substrate, the level of accuracy reduces by roughly two orders of magnitude. 

Figure 6.3 illustrates the magnitude of the magnetic vector potential versus radial distance 
over the range of 10~ 8 < p < 1 for substrate thickness of h = 10" 3 and 10” 1 as computed from 

(6.5). To demonstrate the efficiency of the numerical routine, the total number of function 
evaluations (i.e., the number of evaluations of the integrand) required to compute each point in 
Fig. 6.3 is illustrated in Fig. 6.4. This includes the contour path integration and the method of 
averages. It is seen that as the radial distance decreases, the number of points increases. As seen 
in Fig. 6.2, for small values of p, small argument approximations can be used to accurately 

evaluate the Sommerfeld integral. Nevertheless, for the error criterion specified, at most 1800 
evaluations was required. Thus, the computation of Fig. 6.3 can be done in seconds. 

The same example is repeated for the electric scalar potential (6.6). These results are again 
illustrated in Figs. 6.4 through 6.6, with similar results as the magnetic vector potential. 
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Fig. 6.2. Percent error between magnetic vector potential and a small argument approximation. 



Fig. 6.3 Magnitude of magnetic vector potential. 


DAAG55-98-1-0460 


32 


The University of Kentucky 






STIR: High-Order Sparse Matrix Nystrom Algorithm to the Analysis ofPlanar_ Microwave Circuits - Final Re P ort 



Fig. 6.4 Number of function evaluations for computing potentials in Fig. 6.3. 



Fig. 6.5 Percent error between electric scalar potential and the small argument approximation. 
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Fig. 6.6 Magnitude of magnetic vector potential. 



Fig. 6.7 Number of function evaluations for computing potentials in Fig. 6.6. 
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6.3 Interpolation of the Green's function 

Since the proposed numerical scheme is a high-order solution of the MPIE, the Green's 
function must be evaluated to a controllable accuracy. In the previous section, a means to 
perform the Sommerfeld integrals to controllable accuracy was presented. Unfortunately, it 
would be computationally prohibitive to perform the Sommerfeld integration for each value of 
p required by adaptive quadrature evaluation for the impedance matrix. Rather, it is much more 

efficient to precompute the magnetic vector and scalar potential at discrete points a priori and 
then compute the Green's function for arbitrary values of p via an interpolation scheme. It must 

be realized that to maintain high-order accuracy, the interpolation scheme must be able to 
interpolate the Green's function as well as its derivative to controllable accuracy. This is the 
topic of this section. 

Due to the azimuthal symmetry of the Green's function, the Sommerfeld integrals can be 
tabulated as a one-dimensional array if all the conductors are horizontal and are located at the 
same vertical height. Otherwise, a two-dimensional array need be computed and stored. The 
discrete distances are sampled based on a logarithmic scaling along the radial direction over a 
pertinent range to resolve the large scale behavior of the Green's function. The Green's function 
is thus interpolated based on this logarithmic spacing of points. 

It has been found that polynomial interpolation of the Green's function, e.g., Lagrange 
interpolation polynomials, can not provide controllable accuracy without a very dense sampling 
of points. This is undesirable since it adds to the precomputation cost, as well as the cost of each 
interpolation. A more efficient means of interpolation is a rational function interpolation [11]. 
Specifically, it has been found that a rational function interpolation can estimate the number of 
data points are necessary to controllable accuracy. Furthermore, efficient schemes of evaluating 
the rational function interpolation in multiple dimensions have been developed [15]. 

To this end, we have applied a multidimensional rational polynomial interpolation of the 
layered Green's function. The interpolation scheme is efficiently applied using the Bulirsch- 
Stoer algorithm [16], which is a recursive scheme. Again, the data is logarithmically scaled. 
The Green's function is then computed via interpolation by using local reference points. 
Specifically, for a given value of p , the discrete interpolation points used are determined by a 

fixed window centered about p . 

As an example, the magnetic vector potential in Fig. 6.4 with substrate thickness 10“'* is 
interpolated via a rational function interpolation. To this end, the magnetic vector potential was 
computed at 222 logarithmically spaced points over the range of p e [lO“ 7 ,l] X. Subsequently, 

the magnetic vector potential was then computed at logarithmically spaced points over this same 
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range not coincident with the original data points. The magnetic vector potential was then 
predicted using both rational function interpolation and using the Sommerfeld integration. The 
relative error of the interpolated values versus the integration is illustrated in Figs. 6.8 and 6.9 for 
the real and the imaginary parts of the vector potential, respectively. The interpolations were 
done using both a 5 point and a 10 point window. It is observed that for this range of values, the 
5 point interpolation provides 2 digits or more of accuracy, and the 10 point window provides at 
least 6 digits of accuracy relative to the original integral. It is also observed that the error is more 
accurately predicted for small values of p, while larger errors are encountered for larger values 
of p . This is mainly due to the fact that for larger values of p the magnetic vector potential is 
more oscillatory. As a result, denser sampling denser becomes necessary for larger values of p . 
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Fig. 6.8 Relative error in the real part of the magnetic vector potential predicted by 5 
point and 10 point rational function interpolation over the range of 

[10 ~ 7 A,U]. 



Fig. 6.9 Relative error in the imaginary part of the magnetic vector potential 

predicted by 5 point and 10 point rational function interpolation over the 

range of [l0~ 7 /i, l/l J. 
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7. Conclusions 


In this report, we have summarized the 6 month progress funded under this STIR towards a 
high-order fast solution of the MPIE for the analysis of microwave circuits printed in a layered 
media. To this end, a locally corrected Nystrom solution of the electric field integral equation 
(EFIE) and the magnetic field integral equation (MFIE) for the three-dimensional ffee-space 
Green's function has been developed. Through numerical examples, it was demonstrated that the 
LCN technique is high-order convergent. Specifically, exponential convergence is realized by 
simply increasing the quadrature order of the underlying Nystrom scheme. Furthermore, it was 
demonstrated that high-order convergence is realizable with only modest increases in 
computational time. It was also demonstrated that high-order convergence can be realized for 
metallic structures with sharp edges. This was achieved through the use of a Gauss-Jacobi 
quadrature rule. For flat edges, a Gauss-Jacobi quadrature rule exactly integrates the edge 
singularity and will provide exponential convergence. For curved edges, the singularity is not 
represented exactly. However, it is closely approximated, and the LCN method will provide 
rapid convergence to engineering accuracy. 

The advantages of a high-order scheme such as the locally-corrected Nystrom scheme is that 
solutions to engineering accuracy (1 %) can be realized with less resources than a traditional 
low-order method. More importantly, the LCN method can also offer an error estimate of the 
final solution with little additional computational overhead. 

The second important thrust area presented is in the efficient computation of the layered 
Green's function. To this end, a contour path originally suggested by Gay-Balmaz and Mosig 
[13] is utilized along with the method averages [12, 14] to efficiently compute the Sommerfeld 
integrals to controllable precision. It was then shown that a rational function interpolation of the 
data then provides a highly efficient means of interpolating the data to controllable accuracy in 
one and two dimensions. 
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